% Setting up the parameters.

N1=[22,83,119,12,2];
N2=[8,58,143,23,3];
N3=[10,36,153,34,9];
N4=[8,42,137,45,11];
N5=[4,28,136,44,18];

total = [0;0;0;0;0];
total(1) = sum(N1);
total(2) = sum(N2);
total(3) = sum(N3);
total(4) = sum(N4);
total(5) = sum(N5);

pObs = [N1/total(1);N2/total(2);N3/total(3);N4/total(4);N5/total(5)]
num2str(100*pObs,3)

% ML estimate

% [x,y] = fminsearch(@logL,[.3,0.4,.1])

% pML = pmat(x);

% I'm using the matrix that Aldo sent
%A = load('ML_Electricity.txt')
A = dlmread('ML_Electricity_full.txt');
pML = A
num2str(100*pML,3)



dist_Obs = sum(sqrt((((pObs - pML).*(pObs - pML))')))/5

% This creates the distance measure for which we want to derive the
% distribution
dist_Obs_overall = sum(sum(sqrt((pObs - pML).*(pObs - pML))))/25
num2str(100*dist_Obs_overall,3)

% Create simulations for each row

addup =[1,1,1,1,1;0,1,1,1,1;0,0,1,1,1;0,0,0,1,1;0,0,0,0,1];

% This creates the upper and lower bounds for the quintiles in each row

upper = pML*addup;
lower = [zeros(5,1),upper(:,1:4)];

% This sets the number of simulations

Sims = 100000

% Initialize the variables
dist_Sim_overall = zeros(1,Sims);
dist_Sim_1 = zeros(5,Sims);

for i = 1:Sims


N_Sim = zeros(5,5);

% The following two loops create one experiment
for j = 1:5

    %total(j) is the number of subjects in a row
for k = 1:total(j) 
    m = rand(1);
    
    % Here is what this line does: m is uniform between 0 and 1. upper(j,:) and
    % lower(j,:) contain the cutoffs for each quintile. The code adds one
    % observation to the quintile in which m falls.
    N_Sim(j,:) = N_Sim(j,:) +(lower(j,:)<m & m <=upper(j,:));
end;
end;

% This converts everything into probabilities
pSim = [N_Sim(1,:)/total(1);N_Sim(2,:)/total(2);N_Sim(3,:)/total(3);N_Sim(4,:)/total(4);N_Sim(5,:)/total(5)];

% Sum up the distances from the simulation and store it in a huge matrix.

dist_Sim_overall(1,i) = sum(sum(sqrt((pSim - pML).*(pSim - pML))))/25;
dist_Sim(:,i) = sum(sqrt((pSim - pML).*(pSim - pML))')/5;

end;


% This generates the critical values for the simulated distribution for
% the following p-values: p = 0.2, 0.1, 0.05, 0.01, 0.005, 0.001

quantile(dist_Sim_overall,[.001 .005 .01 .1 .2 .3 .4 .5 .6 .7 .8 .9 .95 .99 .995 .999])

pvalue=sum(dist_Sim_overall>dist_Obs_overall)/Sims
